Load in Precip data

load('data/DailyP.RData')

head(daily_p)
## # A tibble: 6 x 8
##   date       station   lon   lat mean_temp min_temp max_temp daily_p
##   <date>     <chr>   <dbl> <dbl>     <dbl>    <dbl>    <dbl>   <dbl>
## 1 2018-10-01 04V     -106.  38.1      55.7     46.4     71.6    0   
## 2 2018-10-01 0CO     -106.  39.8      38.5     33.8     46.4    0.04
## 3 2018-10-01 20V     -106.  40.1      54.7     39.2     69.8    0   
## 4 2018-10-01 2V5     -102.  40.1      50.9     44.6     68      0   
## 5 2018-10-01 2V6     -103.  40.1      48.7     42.8     68      0   
## 6 2018-10-01 33V     -106.  40.8      54.4     37.4     66.2    0

Get Elevation Data

unique_asos <- daily_p %>%
  distinct(lon, lat, station)  %>%
  st_as_sf(., coords = c('lon','lat'), crs = 4326) %>%
  get_elev_point(.)

Get Monthly P Averages

unique_asos <- st_read('data/unique_asos_elev.gpkg')
## Reading layer `unique_asos_elev' from data source 
##   `C:\Users\kafre\Documents\CSU_R_Projects\ESS580\Intro to Water Data\metals_interp\data\unique_asos_elev.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 134 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -108.7593 ymin: 37.1515 xmax: -102.241 ymax: 40.7503
## Geodetic CRS:  WGS 84
monthly_p <- daily_p %>%
  mutate(month = month(date)) %>%
  group_by(month, station) %>%
  summarize(monthly_p = sum(daily_p)) %>%
  left_join(unique_asos) # grab elevation data
## `summarise()` has grouped output by 'month'. You can override using the `.groups` argument.
## Joining, by = "station"

Look at monthly P

ggplot(monthly_p, aes(x = elevation, y = monthly_p, color = month)) + 
  scale_color_viridis_c() + 
  geom_point()

Getting Monthly Means of means, mins, maxes.

monthly_t <- daily_p %>%
  mutate(month = month(date)) %>%
  group_by(month, station) %>%
  dplyr::select(-lon,-lat) %>%
  summarize(across(where(is.numeric), mean, na.rm = T)) %>%
  left_join(unique_asos,.) 
## `summarise()` has grouped output by 'month'. You can override using the `.groups` argument.
## Joining, by = "station"

Temp vs Elevation

ggplot(monthly_t, aes(y = mean_temp, x = elevation, color = month)) + 
  geom_point() + 
  scale_color_viridis_c()

Pick a month (summer months are safer)

#filter data for June

Build IDW precip or elevation for state for that month

#Using temperature for June
unique_asos_2163 <- st_transform(unique_asos,crs = 2163)

#convert to stars object
co_box <- st_bbox(unique_asos_2163) %>%
  st_as_stars(dx = 1000)

#filter for June
june_t <- monthly_t %>%
  dplyr::filter(month == 6) %>%
  st_transform(., st_crs(co_box)) %>%
  na.omit(.)

#conduct inverse distance weighting
interp_basic = idw(mean_temp~1, june_t, co_box) %>%
  dplyr::select(1)
## [inverse distance weighted interpolation]
#create map
tm_shape(interp_basic) + 
  tm_raster(palette = 'Reds', style = 'cont')

Plot this data

ggplot(june_t, aes(x = elevation, y = mean_temp)) + 
  scale_color_viridis_c() + 
  geom_point() +
  ylab('Temperatureon, June 2019') +
  xlab('Elevation (Meters)')

Build IDW with elevation for state for that month including elevation as a predictor

Hint! Use get_elev_raster

#using the 'raster' library
library(raster)

#get elevation raster
ras <- get_elev_raster(unique_asos, z = 7) %>%
  raster::crop(.,unique_asos)

#convert to STARS object
co_stars <- st_as_stars(ras)

#change row name to elevation
names(co_stars) <- 'elevation'

#filter temperature to june
june_t <- monthly_t %>%
  filter(month == 6)


#make a plot
ggplot(june_t, aes(elevation, mean_temp)) + 
  geom_point() + 
  geom_smooth(method = 'lm')

#conduct interpolation
interp = gstat::idw(mean_temp~elevation, 
                    june_t, 
                    co_stars) %>%
  dplyr::select(1)
## [ordinary or weighted least squares prediction]

Make a Map of that

tm_shape(interp) + 
  tm_raster(palette = 'plasma', style = 'cont')

mapview(interp)

Compare both maps to PRISM approach for your month

How close do our simple approaches come to reproducing prism maps at https://www.prism.oregonstate.edu/recent/monthly.php ?

Ultimately, the PRISM map is fairly close to our generated map, with a different color scheme.